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computations  needed  for  prediction,  design  and/or  analysis  applications.  In  all  cases  for  which  an  associated  flow 
rule  is  deemed  inadequate  to  fit  observed  material  behavior,  an  approach  other  than  rate-independent  plasticity 
must  be  sought  in  order  to  represent  such  behavior  in  a  reliable,  self-consistent  and  rational  manner. 
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SUMMARY 


Many  ground  shock  calculations  performed  for  the  Defense  Nuclear  Agency  utilize  rate-in¬ 
dependent  plasticity  models  with  non-associated  (i.e.,  non-normal)  flow  rules  to  represent  the 
behavior  of  geological  materials.  In  such  models  the  yield  surface  is  used  to  limit  the  shear  and/or 
tensile  stresses,  while  the  flow  rule  is  used  to  independently  define  the  direction  of  plastic  straining 
(thereby  limiting  the  dilatancy).  The  behavior  of  constitutive  formulations  of  this  kind  in  dynamic 
applications  is  analyzed  in  this  report.  By  means  of  a  “constructive”  proof,  it  is  demonstrated 
mathematically  that  these  models  lead  to  multiple,  and  therefore  spurious,  solutions  in  wave 
propagation  situations;  in  other  words,  uniqueness  of  solution  breaks  down  fo1'  dynamic  in¬ 
itial/boundary  value  problems.  Although  a  simple  class  of  multiple  solutions  is  presented  to 
complete  the  proof  for  a  very  simple  type  of  problem,  the  conclusion  is  quite  general;  it  is  relevant 
to  all  of  the  practically  used  ground  shock  (and  some  of  the  structural)  plasticity  models.  The  result 
implies  that  normality  of  flow  is  a  necessary  (as  well  as  a  sufficient)  condition  for  uniqueness  in 
dynamic  applications  of  any  rate-independent  plasticity  model. 

Obviously,  these  findings  also  imply  an  inherent  lack  of  robustness  in  numerical  analyses 
based  on  non-associated  plasticity.  Because  the  use  of  these  models  is  widespread  in  ground  shock 
(and  structural)  calculations,  the  fact  that  they  can  produce  multiple  (and  therefore  spurious) 
solutions  is  significant;  it  casts  serious  doubt  on  the  validity  of  the  specific  results  of  any  (and  every) 
calculation  based  on  this  kind  of  constitutive  representation.  Therefore  it  is  strongly  recommended 
that  ground  shock  calculators  promptly  abandon  their  reliance  on  such  models  in  all  work.  To 
describe  the  situation  most  bluntly,  these  models  do  not  provide  a  rational  basis  for  the  computations 
needed  for  prediction,  design  and/or  analysis  applications.  In  all  cases  for  which  an  associated  flow 
rule  is  deemed  inadequate  to  fit  observed  material  behavior,  an  approach  other  than  rate-independent 
plasticity  must  be  sought  in  order  to  represent  such  behavior  in  a  reliable,  self-consistent  and  rational 
manner. 

In  view  of  the  fact  that  non-associated  plasticity  models  have  already  been  extensively  used 
in  analyses  for  DNA,  it  is  recommended  that  some  effort  be  made  to  quantitatively  assess  the 
consequences  of  their  use.  In  particular,  the  extent  to  which  the  model  assumptions  and  parameters 
can  affect  the  numerical  solutions  to  various  types  of  problems  should  be  investigated. 

Ai?C; 

tT’v 

#nc  QI3  AL;  •  I A  ox-  iiu  i  Au  & 

V  . 

iii 


CONVERSION  TABLE 


CONVERSION  FACTORS  FOR  U.S.  CUSTOMARY  TO 
METRIC  (SI)  UNITS  OF  MEASUREMENT 


MULTIPLY  . *~BY  . TO  GET 

TO  GET  . BY  . DIVIDE 


bar 

1.000  000XE  +  2 

kilopascal  (kPa) 

degree  Farenheit 

(1^-32)/!^ 

degrees  Celsius 

foot 

3.048  000  X  E  l 

meter  (m) 

foot-pound-force 

1.355  818 

joule  (J) 

inch 

2.540  QOOXE-2 

meter  (m) 

kip  (1000  Ibf) 

4.448  222 XE  + 3 

newton  (N) 

kip/inch  2  (ksi) 

6.894  757  XE  + 3 

kilopascal  (kPa) 

kips  per  foot 

14.5932 

kilonewtons  per  meter 

kips  per  sq.  ft.  (ksf) 

4.788  X  E  +  3 

pascal  (Pa) 

ksi 

6.894  757XE  +  6 

pascal  (Pa) 

pound 

0.4536 

kilogram  (kg) 

pound-force  (lbs  avoirdupois) 

4.448  222 

newton  (N) 

pound-force  inch 

1.129  848  X  £  - 1 

newton-meter  (N  •  m) 

pound-force/inch 

1.751  268XE  +  2 

newton/meter  (N/m) 

pound-force/foot 2 

4.788  026  XE- 2 

kilopascal  (kPa) 

pound-force/inch 2  (psi) 

6  894  757 

kilopascal  (kPa) 

pressure  (psi) 

6.894  757  XE  4*3 

pascal (Pa) 

square  inch 

6.452  XE- 4 

meter 2  (m  2 ) 

square  foot 

9290  X  E  -  2 

meter 2  (m 2 ) 

slug 

1.459  390XE  +  1 

kilogram  (kg) 

TABLE  OF  CONTENTS 


Section  Page 

i 

SUMMARY .  iii 

CONVERSION  TABLE .  iv 

1  INTRODUCTION . I 

2  RATE-INDEPENDENTNON-ASSOCIATED  PLASTICITY .  2 

3  DERIVATION  OF  THE  TANGENT  STIFFNESS  TENSOR .  5 

4  WORK  DONE  IN  AN  INFINITESIMAL  LOAD-UNLOAD  CYCLE  7 

5  ELASTIC-PLASTIC  WAVE  PROPAGATION .  11 

6  FAST  PLASTIC  WAVES  AND  NON-UNIQUENESS .  20 

7  CONCLUSION .  27 

APPENDIX 

DERIVATION  OF  THE  STRAIN  INCREMENTS  SATISFYING 
INEQUALITY  (4.7)  .  29 


v 


SECTION  1 
INTRODUCTION 


To  represent  geological  materials  in  ground  shock  calculations,  analysts  often  use  rate-in¬ 
dependent  plasticity  models  which  are  non-associated,  i.e.,  which  utilize  plastic  strain  increments 
not  directed  along  the  normal  to  the  yield  surface.  Although  such  models  are  often  fit  to  quasistatic 
material  tests,  they  are  used  in  dynamic  ground  shock  analyses  in  a  rate-independent  formulation, 
with  only  the  values  of  the  material  constants  being  adjusted  to  account  for  dynamic  effects.  In  such 
models  the  yield  surface  is  used  to  limit  the  shear  and  tensile  stresses,  while  the  flow  rule  is  used  to 
independently  define  the  direction  of  plastic  strain  (and  thereby  limit  the  dilatancy  exhibited  by  the 
model). 


The  wave  propagation  behavior  of  these  models  is  analyzed  in  this  report.  Although  a 
considerable  amount  of  test  data  is  consistent  with  non-associated  plastic  flow  (Lade  1988)1,  it  is 
shown  below  that  the  assumption  that  this  applies  for  all  strain  rates  is  not  appropriate  in  dynamic 
analysis.  Viscosity,  size  effects,  etc.,  in  real  materials  may  affect  the  flow  direction  for  strain  rates 
of  importance;  in  fact  such  effects  may  preclude  any  plasticity-based  formulation  from  working 
well.  The  present  analysis  proves  that  the  rate-independent,  non-associated  plasticity  used  in  practice 
is  not  rational  (i.e.,  it  leads  to  improperly  posed  problems  possessing  multiple  -  rnd  therefore 
spurious  -  solutions)  whenever  inertia  effects  are  included  (as  is  always  the  case  for  wave 
propagation  applications). 

This  report  presents  the  general  equations  of  rate-independent  plasticity  in  Section  2.  This 
is  followed  by  a  derivation  of  the  plastic  tangent  stiffness  tensor  in  Section  3,  an  analysis  of  certain 
infinitesimal  load-unload  cycles  in  Section  4,  and  a  derivation  of  the  general  equations  governing 
elastic-plastic  wave  propagation  in  Section  5.  Finally,  the  results  of  Sections  2-5  are  utilized  in 
Section  6  to  construct  a  simple  initial/boundary  value  problem  f  ,t  which  the  existence  of  multiple 
solutions  is  easily  demonstrated.  It  has  long  been  known  that  associated  flow  is  a  sufficient  condition 
for  uniqueness  in  plasticity.  This  report  proves  that  it  is  also  a  necessary  condition  for  uniqueness 
for  the  rate-independent  models  used  in  practical  analyses  of  ground  shock. 


1  P.V.  Lade,  “Effects  of  Voids  and  Volume  Changes  on  the  Behavior  of  Frictional  Materials”, 
Inti.  J.  Num.  Anal.  Meth.  Geomech.,  Vol.  12,  No.  4,  p.351-370,  1988. 


SECTION  2 

RATE-INDEPENDENT  NON-ASSOCIATED  PLASTICITY 


Non-uniqueness  will  be  proven  by  demonstrating  that  multiple  solutions  can  occur  during 
an  infinitesimally  short  time  interval  dt.  During  such  an  interval  the  changes  in  the  displacement 
field  itself  are  infinitesimal.  Therefore  we  may  linearize  the  geometric  treatment  without  loss  of 
generality,  provided  that  we  adopt  a  reference  configuration  corresponding  to  the  (possibly 
deformed)  current  geometry  at  the  instant  of  interest.  If  we  adopt  such  an  “updated  Lagrangian” 
approach,  then  we  may  utilize  the  linearized  “small  displacement”  kinematic  strain  decomposition 
during  the  interval  dt. 


■E  P 
E  •  =  £  -  f  £ 
c  IJ  ^lj  r 


(2.1) 


where  the  components  of  the  strain  rate  tensor,  d  noted  by  are  measured  with  respect  to  a 
reference  configuration  which  represents  the  geometry  at  the  instant  of  interest.  The  superscripts 
denote  the  elastic  and  plastic  strain  rates,  and  the  elastic  strain  is  related  to  the  stress  through  the 
elastic  stiffness.  (In  this  report  we  adopt  the  engineering  mechanics  convention  that  stresses  and 
strains  are  positive  in  tension). 


Aside  from  the  kinematic  decomposition  of  the  deformation  into  elastic  and  plastic  strains, 
all  rate-independent  plasticity  models  are  characterized  by  means  of  four  main  features.  These  are 
the  elastic  behavior  (to  be  discussed  later),  the  yield  condition,  hardening  (if  any),  and  the  flow  rule. 
The  yield  condition  can  be  written  most  generally  as  f(G;j  ,  K*)=0,  where  are  the  components  of 
the  stress  tensor  and  the  Kk,k=l  ...n,  denote  the  n  hardening  parameters  in  the  model  (which  are 
absent  in  the  special  case  of  ideal  plasticity).  The  yield  condition  may  be  considered  a  hypersurface, 
called  the  yield  surface,  in  a  six-dimensional  stress  space.  (Each  point  in  such  a  space  corresponds 
to  a  stress  state,  and  the  yield  surface  contains  those  points  corresponding  to  stresses  at  which  plastic 
deformation  can  occur).  Without  loss  of  generality  we  can  associate  values  /<0  with  elastic  states 
“within”  the  yield  surface,  while  the  region  /> 0  “outside”  the  yield  surface  is  excluded  by  the  theory . 
This  exclusion  implies  that  /= 0  for  as  long  as  plastic  yielding  occurs,  so  that  differentiation  of  the 
yield  condition  with  time  gives  /  =  0,  the  “plastic  compatibility  condition” 


JL*„  +  iL,.0 

do,  ‘J  c)K*  * 


(2.2) 
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which  must  be  satisfied  during  plastic  response.  (In  this  report  we  adopt  the  notation  =  —  and 

utilize  the  Einstein  summation  convention  concerning  subscripted  indices,  except  where  explicitly 
indicated  otherwise). 

Because  hardening  (if  any)  proceeds  according  to  specified  “hardening  rules”  as  plastic 
straining  occurs, 

(23) 

where  the  ly  are  coefficients  used  to  represent  the  hardening  behavior.  Equation  (2.2)  therefore 
becomes 


nIj  °ij  =  Hij  ej 


(2.4) 


in  which 


jV 

9ki . 


(2.5) 


It  should  be  noted  that  the  tensor  nfj  corresponds  to  the  outward-drawn  “normal”  to  the  yield  surface 

Y  y 

in  stress  space  and  possesses  the  symmetry  n{j  =  n /l-  (because  of  the  symmetry  of  Cy). 

To  complete  the  me  del,  the  flow  rule  must  be  specified.  This  is  done  in  terms  of  a  “direction” 
of  plastic  strain  ny  so  that  the  flow  rule  can  be  written  generally  as 


(2.6) 


in  which,  without  loss  of  generality,  is  chosen  so  that  k> 0  during  plastic  deformation  (k= 0  during 

elastic  deformation),  and  the  magnitudes  of  n y  and  X  are  scaled  so  that 


P  P  Y  Y 
nijnij  nijnij 


(2.7) 
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P  Y 

If  the  "normality  condition"  n i}  =  is  satisfied,  the  flow  rule  is  said  to  be  normal,  or  associated; 

P  Y  P  P 

otherwise  n y  nl}  and  the  flow  rule  is  said  to  be  non-associated.  It  should  be  noted  that  ntJ  =  n  t 

(because  of  the  symmetry  of  ey). 

Summarizing,  rate-independent  plasticity  models  postulate  that  material  deformation  con¬ 
sists  of  elastic  and  plastic  components,  and  all  such  models  exhibit  four  features: 

a)  an  elastic  stiffness  tensor  that  defines  the  functional  relationship  between  the  stress 
and  the  elastic  strain, 

b)  a  yield  function  that  defines  the  plastic  states  of  the  material,  and  which  is  defined 
in  terms  of  the  stress  as  well  as  hardening  parameters  (except  in  the  case  of  ideal 
plasticity), 

c)  hardening  rules  which  specify  the  way  the  hardening  parameters  evolve, 

d)  a  flow  rule  that  specifies  the  “direction”  of  plastic  straining. 

Particular  choices  for  the  arbitrary  functions  in  a-d  define  the  individual  models  which  are  used  in 
practical  applications,  and  which  make  up  the  broad  class  of  models  examined  in  this  report. 
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SECTION  3 

DERIVATION  OF  THE  TANGENT  STIFFNESS  TENSOR 


We  begin  our  analysis  of  the  wave  propagation  characteristics  of  plasticity  models  by 
deriving  the  general  form  of  the  tangent  (elastic-plastic)  stiffness  tensor.  The  incremental  relation¬ 
ship  between  the  stress  and  strain  rates  for  elastic  response  is 

®ij  ~  Cijkl  ^ kl  (3.1) 

where  C,yW  are  the  components  of  the  (fourth  rank)  elastic  stiffness  tensor.  It  is  well  known  'rat 
CyM  possesses  the  symmetries  C,^/  =  Cy#  =  Cj =  C*/y.  In  the  plastic  case  substituting  eq.  (2. 1 ) 


gives 


®mn  Cmnpq  i^pq~^pq) 


Multiplying  this  equation  by  n ^  and  substituting  eqs.  (2.4)  and  (2.6)  leads  to 

mnnmn  ~  nmrS'mnpq^'pq~^npq)  ~  ^~^pqnpq 

where  some  dummy  indices  have  been  interchanged.  Equation  (3.3)  can  be  solved  for  X, 


yr  * 

^  =:  Qnmn^mnpq  ^ pq 


in  which 


C mnpqKpq  +  H pqnpqj 

In  standard  dynamic  computational  procedures  the  constitutive  behavior  is  utilized  to 
calculate  the  stresses  from  strains  (which,  in  turn,  are  computed  from  the  field  equations).  We 
therefore  consider  E,y  to  be  the  input  here  and  6y  to  be  the  desired  output  quantity.  According  to  eq. 
(3.1)  the  quantity  Cmnpq  epq  represents  the  rate  of  change  of  stress  in  the  case  of  elastic  response 

(i.e.,  epq  -  0),  so  the  first  of  eqs.  (2.5)  implies  that  the  scalar  nmn  Cmnpq  tpq  represents  /  ,  the 
hypothetical  rate  of  change  of  the  value  of  the  yield  function  in  that  case.  Because  the  definition  of 

•  c 

/  implies  that  a  negative  value  of  /  would  indeed  be  consistent  with  elastic  behavior,  plasticity 
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*  F 

can  occur  only  for  non-negative  values  of  /  (otherwise  the  stress  response  would  not  be  unique). 
Furthermore,  if  /l  eq.  (3.4)  indicates  that  X  =  0,  so  that  iy  =0  from  eq.  (2.6).  Therefore,  /E=0 
also  implies  that  the  response  is  elastic  (this  situation  is  called  neutral  loading).  For  these  reasons, 
plastic  yielding  (£-,-  *  0)  will  occur  if  and  only  if 

nY  C  £  >0  (3.6) 

nmn  '•'mnpq  cpq  w 


P 

Because  X>0  in  this  case,  eq.  (3.4)  implies  that  Q  >0.  For  consistency,  therefore,  npq  must  be  defined 
so  that  the  denominator  of  eq.  (3.5)  is  always  positive.  Substituting  eqs.  (2. 1 ),  (2.6)  and  (3.4)  into 
(3.1)  gives  the  incremental  stress-strain  relation 


P  _Y 


®ij  ~  Cijkl  &kl  QCijkl  nkl  nmn^~mnpq  &pq 


(3.7) 


which  (after  interchanging  the  dummy  indices  k,l  and  p,q  in  the  last  term)  gives 

•  .  p  y 

®ij  =  Cijkl  ^ kl  ~  Q^ijpq  npq  nmn^mnkl  E kl 

Invoking  the  symmetry  Cmnkl  -  Ckimn  leads  to 


Cjdmn  nmr^j-kl 


(3.9) 


The  quantity  in  the  parenthesis  of  eq.  (3.9)  is  the  desired  tangent  stiffness  modulus,  Ly^.  It 
may  be  written  more  conveniently  as 


p  y 

Lijk\  ~  Cijkl  ~  Qmijmkl 


(3.10) 


where 


mij  —  C ijpq  npq  >  mij  35  ^ijpq  npq  (3. 11) 

Note  that  for  associated  flow  my  =  nty  so  that,  by  virtue  of  eq.  (3. 10),  the  tangent  stiffness 
tensor  possesses  the  “stress-strain”  symmetry  Ly =  Lkly.  For  general  non-associated  flow,  how¬ 
ever,  eq.  (3. 10)  implies  that  the  tangent  stiffness  tensor  does  not  possess  this  particular  symmetry, 
although,  of  course,  Lyk[  =  Ljiki  —  Lylk  in  all  cases. 
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SECTION  4 

WORK  DONE  IN  AN  INFINITESIMAL  LOAD-UNLOAD  CYCLE 


Consider  the  response  of  a  rate-independent  plastic  material  subjected  to  an  infinitesimal 
load-unload  strain  cycle.  (In  order  to  do  this,  we  refine  an  analysis  previously  described  in  (Sandler 
and  Rubin  1987)  ,  which  in  turn  was  partially  based  on  arguments  first  presented  in  (Il'iushin 

a 

1961)  ).  Refer  to  Figure  1,  which  shows  a  two-dimensional  schematic  representation  of  an 
infinitesimally  small  region  of  the  six-dimensional  stress  space.  Point  O  corresponds  to  any  initial 

stress  state  a fj  on  the  yield  surface,  represented  locally  by  the  line  OY,  while  the  inclined  line  OS 

p 

represents  (locally)  the  surface  orthogonal  to  the  “flow  direction”  rty  (this  surface  is  called  the  plastic 
potential  surface).  Let  Aew  be  an  infinitesimal  strain  increment  from  the  state  O  that  causes  the 
response  to  be  plastic.  The  resulting  stress,  shown  in  Figure  1  at  P,  is  given  by 

°ij  =  0%  +  Lijklte-ki  (4- 1  ^ 

/ 

(The  figure  indicates  some  hardening  of  the  yield  surface  during  the  loading  increment,  but  this 
detail  is  not  significant;  for  ideal  plasticity  P  lies  on  OY).  The  unloading  segment  PQ,  which 
corresponds  to  the  strain  increment  -A £*/,  is  elastic  and  results  in  the  final  stress 

cQ  =  G?j-Cijk[&£kl  (4'2) 

The  work  (energy)  required  to  deform  a  unit  volume  of  the  material  through  the  load-unload 
strain  cycle  C  (A followed  by  -Ae^)  is  given  by 


For  the  infinitesimal  cycle  ±A£y  the  Mean  Value  Theorem  can  be  applied  to  obtain 

2  I.  Sandler  and  D.  Rubin,  “The  Consequences  of  Non- Associated  Plasticity  in  Dynamic 
Problems”,  in  Constitutive  Laws  for  Engineering  Materials  -  Theory  and  Applications,  Ed.  by 
C.S.  Desai,  et  al,  Elsevier,  p.  345-352,  1987. 

3  A.A.  Il’iushin,  “On  the  Postulate  of  Plasticity”,  PMM,  Vol.  25,  No.  3,  p.  503-507,  1961. 
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plasticity  model  ,o  n°n““sot,aKd 
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AW  «  +  cfj  ^Eij  +  %Oij  +  o$  )-AEy)  =  -  a f  ^ 


(4.4) 


Referring  back  to  eqs.  (3.10)  and  (4. 1,4.2)  this  becomes 


AW  =  ^Qmfjmli&Eyte.ij 


(4.5) 


Substitution  of  eqs.  (3. 1 1)  finally  gives 


AW  =  i 


=  j  Q^ijp^pqteijJCldrs*  rs&Z-ki'j 


(4.6) 


p  y 

Note  that  for  associated  flow  =  n L  so  that  AW  can  never  be  negative  because  the  two  expressions 

enclosed  in  parenthesis  in  eq.  (4.6)  are  equal  and  Q  is  positive. 

In  the  non-associated  case,  it  is  easy  to  find  strain  increments  Aew  which  satisfy  the  pair  of 
inequalities 

njfiijklteki  >  0  >  nfjCijk!AEki  (4.7) 

For  example,  consider  the  strain  increment  defined  by  the  equation 

Cijklte-ki-nfj-nfj  (4.8) 

In  this  case,  using  eq.  (2.7),  we  can  calculate  the  quantities 


Jr  Ac  -JJ  JJ-1/JJ^JJ  ~>J j\-iL 
nij  '-'ijkl^-kl  ~  nij  nij  nij  nij  ~  J nijnij  +nijnij  *-nij  nij  r 


I  nZ-41 


>0 


(4-9) 


nij  Cijkl^-kl  ~  nij  nij  nij  nij  ~  nij  nij  nij  nij  ~  nij  ^ijkl^kl  <  ® 


(4.10) 


so  that  eq.  (4.7)  will  indeed  be  satisfied.  A  derivation  of  the  complete  (six-dimensional)  set  of 
A eki  for  which  eq.  (4.7)  holds  is  given  in  the  Appendix. 
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Consider  the  geometric  interpretation  of  eq.  (4.7).  Its  first  inequality  merely  implies,  from 
eq.  (3.6),  that  the  strain  increment  Aew  causes  plastic  loading,  with  the  quantity  C^/Ae*/  representing 
the  hypothetical  stress  increment  for  elastic  instead  of  plastic  response.  This  hypothetical  “elastic 
stress  increment”  is  shown  in  Figure  1  as  the  dashed  line  OE.  According  to  the  second  inequality 

in  (4.7),  the  elastic  stress  increment  must  form  a  negative  inner  product  with  Therefore  (4.7) 
merely  requires  that  the  line  OE  lie  within  the  “wedge”  YOS.  We  have  just  shown  that  paths  such 
as  OE  exist  for  non-associated  flow,  but  (4.7)  can  never  be  satisfied  in  the  associated  case,  because 
its  left  and  right  sides,  which  have  opposite  signs,  then  become  equal.  (In  other  words  the  “wedge” 
YOS  in  Figure  1  collapses  into  the  single  line  OY  when  the  flow  rule  is  associated). 

Let  us  now  return  to  eq.  (4.6).  Inequalities  (4.7)  imply  that  the  factors  enclosed  in  parentheses 
in  eq.  (4.6)  are  of  opposite  sign,  so  that  the  work  AW  must  be  negative  (because  Q  >  0).  Therefore, 
instead  of  requiring  that  work  be  dissipated  during  all  strain  cycles  involving  plastic  response,  the 
non-associated  material  model  actually  supplies  energy  whenever  the  cycle  ±Ae,y  satisfies  the 
inequalities  (4.7).  In  the  next  two  sections  we  will  show  how  the  possibility  of  such  “energy 
generation”  causes  uniqueness  to  break  down  and  introduces  spurious  signals  into  dynamic 
initial/boundary  value  problems. 
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SECTION  5 

ELASTIC-PLASTIC  WAVE  PROPAGATION 


In  order  to  understand  the  behavior  of  non-associated  plasticity  models  in  dynamic  problems, 
we  will  first  analyze  the  propagation  of  simple  plane  waves.  Consider  plastic  loading  waves  that 
propagate  with  wave  speed  c  into  material  at  rest  in  a  homogeneous  (uniform)  stress/strain  state  O 
at  yield,  and  produce  motion  behind  the  wave  front  in  a  single  spatial  direction  (which  will  generally 
be  different  from  the  direction  of  propagation  of  the  wave)  for  all  of  the  particles.  The  situation  is 
depicted  in  Figure  2,  in  which  the  unit  vectors  a,  and  ft  represent,  respectively,  the  (constant) 
direction  of  wave  propagation  (normal  to  the  wave  front)  and  the  direction  of  particle  motion 
(displacement,  velocity  and  acceleration).  As  shown  in  the  figure,  the  “local”  coordinate  system, 
dxit  has  its  origin  at  an  arbitrary  point  on  the  wave  front  at  some  arbitrarily  chosen  instant.  Let 
dsscdt-  CLj  dxj  denote  the  infinitesimal  distance  of  point  dxj  behind  the  wave  front  at  a  later  time 
dt.  The  particle  displacement  field  g(<&)ft  due  to  the  wave  produces  the  displacement  field 


ui  ( dxj ,  dt)  =  uf  (dxj)  +  g  (ds)  ft  =  uf  (0)  +  -^dxj  +  g  ( ds )  ft 


(5.1) 


where  uf(dxj)  is  the  initial  displacement  field  in  strain  state  ef,  and  g  is  any  continuous,  non-de¬ 
creasing  function  of  ds  behind  the  wave  front  and  zero  ahead  of  it.  Therefore  g  >  0  and  g  >  0.  The 
remaining  relevant  field  quantities  are  derivable  from  ui  by  means  of  the  equations 


dui  . 

vi  =  ^jj  =  8$i  =  cg  ft 


(5.2) 


dt 


=  *P,  =  cV'0, 


(5.3) 


1  du ;  ditj  n  1 


(5.4) 


°y  =  G?j  +  Lijkl  (£W  ~  eW  )  ~  G?j  +  2Lijkl  (~°-l  P*  “  ak  P/  )g  =  ~  Lijkl  ak  P/  g 
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Plane  Wave  at  Time  t 


\ 


Figure  5.1.  Plane  wave  propagating  in  a,  direction  with  particle 
motion  in  ^direction. 


2 


where  v,-  and  a,  are  the  velocity  and  acceleration,  respectively,  is  the  initial  (uniform)  stress,  and 


,  d 


=  — .  (To  obtain  the  above  equations  the  chain  rule, 


4z  =  dS.dL  .  ^ 

dx;  dsdx ,  ’  3,1(1  8r  dsdt  cg 


(5.6) 


as  well  as  the  “strain”  symmetry  were  invoked).  The  equation  of  motion  is 


po,= 


80y 

dxj 


(5.7) 


in  which  p  is  the  density  of  the  material.  Substituting  eqs.  (5. 3, 5. 5,5.6)  into  eq.  (5.7)  gives 


pc  V&  =  Lijkl  ak  ctj  g' 


(5.8) 


behind  the  wave  front.  By  interchanging  dummy  subscripts  in  the  last  equation  we  may  write 

g”(a-k  Likij  «/  Py  *  P^2Pi  j = g"(ak  Likij  ot,  -  pc2Sy  Jj  =  0  (5-9) 

in  which  5y  is  the  Kronecker  delta  (unity  if  i—j,  zero  if  i*f).  For  a  (non-trivial)  wave  solution  the 
“loading  rate”  g"  will  not  be  zero,  so  that  wave  propagation  requires  that 


6^  =  0 


(5.10) 


in  which 


Ay  s  a*.  Lfoij  a;  ;  and  M  =  pc2  (5.11) 

Equations  (5.10,5.11),  which  must  be  satisfied  for  a  physical  wave  to  exist,  form  an 
eigenvalue  problem.  The  three  eigenvalues  M  of  the  3x3  matrix  A(j  represent  the  values  of  the 

quantity  pc  for  which  wave  propagation  can  occur  (i.e.,  for  which  non-zero  eigenvectors  p,  are 
possible).  These  values  determine  the  wavespeeds  c  which  are  characteristic  of  the  material  (at  its 
current  state  O)  for  the  direction  of  propagation  a,-,  while  the  three  eigenvectors  (5t  correspond  to 
the  directions  of  the  particle  motion  in  the  different  types  of  wave  which  may  occur.  We  will  return 
to  this  point  later. 


Equations  (5.1)  through  (5.11)  apply  as  well  to  elastic  wave  propagation  if  the  tangent 
stiffness  Lyu  is  replaced  by  the  elastic  stiffness  C^.  Define 

Aij  =  a*  Cmj  c i[  —  ak  Ckjj\  <X/  —  ctk  Cjikj  o iq  —  (X/  Cjkn  ctk  =  ctk  Cjkn  ctj  =  Ajj  (5.12) 

so  that  the  eigenvalue  problem  corresponding  to  elastic  waves  is 

AfjPj  =  ME  p,-  (5.13) 

where  ME  is  an  eigenvalue  of  AE.  Because  Cyki  is  positive  definite,  the  matrix  AE  is  also  positive 

definite  and  symmetric  in  view  of  eq.  (5.12);  therefore  it  has  three  real  and  positive  eigenvalues, 

F  th 

denoted  by  Mk  and  a  corresponding  set  of  orthonormal  eigenvectors  z*.  If  we  define  zki  as  the  i 

component  of  z*,  we  can  write 

AE  zkj  -  Aff  zki  (no  sum  on  it)  ;  where  Aff  >  Aff  >  A/f  >  0  and  zk(Zu  =  5*/  (5-14) 

Therefore  three  different  elastic  waves  are  possible  (for  k=  1,2  or  3)  in  which  (3,  =  zw.  Also,  note 

£ 

that  -z*  can  be  chosen  instead  of  z*  as  an  eigenvector  of  Ay;  we  shall  choose  each  z*  in  such  a  way 
that  zkflLi  ^  0. 


As  an  example,  let  us  consider  linear  isotropic  elasticity.  In  this  simple  case  the  elastic 
stiffness  is  expressible  as 


2  G 


cikij  ~(K-  o  )5(fc5y  +  G($ifikj  +  &i/>ki) 


(5.15) 


where  K  >  0  and  G  >  0  are  the  bulk  and  shear  moduli  of  the  material,  so  that 


2  G 


~  a k  Cikij  « i  =  (K~  j  +  G  a ,•  oc,  +  G  5 tj  ak  ak 


=  (K  +  ^)aiaj  +  Gbij 


(5.16) 


7  E 

because  ak  ak  =  la*l  =  1 .  Multiplying  eq.  (5. 1 6)  by  zkj  and  subtracting  Mk  zk(  from  both  sides  gives 
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(5.17) 


Ay  Zfy  -  Aff  zy  =  (AT  +  — )  a,  (aj  zkj  )  +  {G  -  Mk)zy  =  0  (no  sum  on  k) 
which,  after  multiplying  by  a,,  leads  to 


Ay  «i  Zkj  ~  ai  Zki  =(K  +  ^y-Mf)(al  zy )  =  0  (no  sum  on  k) 


(5.18) 


E  4(j 

so  that  either  Af*  =  K +  —  or  a(-  z^  =  0  must  be  satisfied.  For  the  first  possibility,  which  -  for  reasons 

that  will  become  apparent  shortly  -  we  will  label  k=l,  substitution  for  A/f  into  eq.  (5.17)  implies 

that  a,  (a.j  z\j )  =  zy  so  that  zj,-  oq  (tty  z^ )  =  (a;  zy )  must  be  equal  toz1(  zy  =  j  =  1.  Because  we 
have  chosen  zy  so  that  z](-a,  >0,  ZyCti  =  +1,  implying  that  zi,  =  at;  such  a  wave  is  said  to  be  a 
longitudinal,  primary  or  P  wave.  For  the  second  possibility,  a,  zy  =  0  implies  that  zy  is  perpendicular 
to  a,-  so  that  such  a  wave  is  a  transverse,  shear,  secondary  or  S  wave;  in  this  case  eq.  (5. 17)  implies 

Af^  =  G.  Because  two  orthogonal  eigenvectors  may  be  chosen  in  the  plane  perpendicular  to  a,, 

Af^  =  G  represents  a  double  root  and  each  line  in  the  plane  represents  one  of  the  possible  directions 
of  polarization  of  the  S  wave.  Therefore,  Z2,  can  lie  in  any  direction  perpendicular  to  a,  with  z3/ 
perpendicular  to  both  a,  and  to  zji-  Summarizing  the  isotropic  linearly  elastic  case,  the  wave  speeds 

must  satisfy  the  inequalities  Aff  =  AT +  —  >  Aff  =  Aff  =  G  >  0,  and  the  direction  of  particle  motion 

Zy  must  coincide  with  a„  while  zn  and  z3l-  can  be  chosen  as  any  two  mutually  perpendicular 

directions  in  the  plane  perpendicular  to  a,.  We  also  note  that  the  strict  inequality  Aff  >  Aff  appears 
to  hold  for  all  elastic  solid  behavior  (anisotropic  as  well  as  isotropic);  we  will  utilize  this  fact  later. 

We  now  return  to  the  analysis  of  plastic  loading  waves,  as  described  by  eqs.  (5.10,5.1 1). 
Applying  eqs.  (3.10)  and  (5.12)  leads  to 


so  that 


Ay  ~  Cikij  o.i -Qak  mik  my  O-i  —  Ay  Q  cq.  my  oq 

Ay  Pj  =  Afj  -  Q  mfk  ak  ( m J  a,  fy)  =  Af  (3, 


(5.19) 


(5.20) 
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Define  a  new  quantity  y,-  as  the  solution  of  the  system  of  linear  equations 


ZpYy*P« 


(5.21) 


so  that,  in  view  of  the  last  of  eqs.  (5.14),  yp  =  hpj  y;  =zpj  zp  yj  =zpi  p,-.  Multiplying  eq.  (5.20)  by  zqi 
and  substituting  eq.  (5.21)  leads  to 


Zqi  A ij  Zpj  Yp  Q  Zgi  mik.  &k  ( mlj  Zpj  Yp  )  —  M  Zqi  Zp{  yp  —  M  5qp  yp 


(5.22) 


From  the  first  of  eqs.  (5.14) 

Zqi  A ij  zpj  —  zqi  Mp  zpi  —  Mp  &pq  (no  sum  on  p) 


(5.23) 


so  that,  eq.  (5.22)  implies 


[  M  Mq  ]  y q  Q  Zqj  mik  ^ Ij  &l  zrj  Y r  ) 


(no  sum  on  q) 


(5.24) 


or 


?4  = 


Q  zQl  mfk  ak 


a 


(no  sum  on  q) 


where  a  represents  the  t;rm  in  parenthesis  in  eq.  (5.24).  Now 

a  s  ~mkl  ak  Zpl  yp  —  ~Cklij  nij  ak  P/ =  ~nij  ^  ^ij^l +  ^ijlkj  ak  P/ 
=  nIj  Cijkl  \  (~<*k  Pi  ~  al  Pit  )  =  nJj  Cijkl  (Ekt  ~  ef/)/ g  >  0 


(5.25) 


(5.26) 


where  the  last  inequality  is  inferred  from  eq.  (3.6)  together  with  g'  >  0  for  the  plastic  loading  wave. 
Furthermore,  from  eqs.  (5.25,5.26), 


3 

—a  =  my  Ctj  Zpj  yp  —  ^(ttj  mg  zpj  )(  zpi  m^  0 ) 
1 


M-Mp 


(5.27) 


Because  a  >  0  for  plastic  loading,  eq.  (5.27)  gives 
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(5.23) 


p\y i  P2 yi  P3 y?>  i  n 

JF  +  F  +  '  E+ 

M- Ml  M-Mf  M-Aff  Q 


in  which 

Pi  =  “  «*  myi  zu  and  yt  =  -ak  mYki  za  (5'29) 

Equation  (5.28)  actually  represents  the  secular  or  characteristic  equation  for  the  eigenvalues 
Mi  in  the  plastic  loading  wave;  it  is  equivalent  to  a  cubic  equation  for  M.  The  quantities  pt  and 
depend  upon  the  direction  of  propagation  and  on  the  elastic  eigenvectors  corresponding  to  that 
direction. 

r 

Let  us  replace  the  unknown  plastic  wave  speed  M  in  eq.  (5.28)  by  x  =  M-M\  and  let 

Mi  -  A/f  =  b  and  Aff  -  A/f  =  c,  where  b  and  c  are  positive.  After  clearing  of  fractions  eq.  (5.28) 
becomes  a  cubic  in  x. 


t3  +  [Q(pi  y,)+h+c]T2  +  [Qpi  yi  (b+c)  +  Qp2  y2c+  QPl  >>3  hlx  +  QP\  >1  bc  =  Q 


(5.30) 


The  sign  of  the  constant  term  in  a  real  cubic  equation  is  always  opposite  to  that  of  at  least 
one  of  the  roots  (which,  of  course,  must  be  real).  Therefore,  if  piyi  is  negative  (which  cannot  occur 
in  the  associated  case,  because  then  pi  =yj),  at  least  one  value  of  x  must  be  positive.  This  m^ans 

that  eq.  (5.28)  possesses  a  real  root  M  >  Aff,  so  that  the  corresponding  plastic  wave  travels  faster 
than  any  of  the  elastic  waves.  It  was  shown  in  (Sandler  and  Rubin  1987)  how  such  a  fast  plastic 
wave  speed  leads  to  non-uniqueness;  that  argument  will  be  reformulated  in  the  next  section,  but 
first  we  need  to  determine  the  conditions  under  which  p\y\  is  negative.  In  order  to  do  this,  note  that 


Comparing  eqs.  (5.31,5.32)  with  eq.  (4.7),  it  is  clear  that  the  condition  p\y j  <  0  is  equivalent  to 
requiring  that  the  “stress  increment”  associated  with  the  “strain  path”  -  ct^n  be  within  the  wedge 
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in  Figure  1.  For  this  reason,  any  value  of  a,  for  which  pi  and  y\  have  opposite  signs  can  be  used  to 
construct  a  "strain  path"  in  the  wedge,  and  represents  a  direction  of  propagation  for  which  plastic 
waves  travel  faster  than  elastic  waves. 

In  order  to  prove  that  such  propagation  directions  ot,  always  exist  for  the  non-associated 
flow  models  used  in  practical  applications,  we  must  introduce  a  new  concept  and  definition.  Note 
that  in  all  such  practical  plasticity  models  the  yield  surface  limits  the  shear  stresses  at  many  (if  not 
all)  stress  states,  while  the  flow  rule  permits  the  relatively  large  shear  strains  associated  with  yield 

or  failure.  Let  us  describe  the  tensor  m ,y  (or  m,y)  as  “shcar-like”  if  at  least  one  of  its  eigenvalues  is 
positive  and  at  least  one  is  negative.  In  other  words,  one  of  the  principal  directions  of  a  shear-like 
tensor  represents  tension  (or  extension)  while  another  represents  compression  (or  contraction).  With 

y 

this  definition  we  can  say  that  all  practical  plasticity  models  have  many  stress  states  for  which 

p 

and/or  m,y  is  shear-like. 

Let  mkl  be  any  shear-like  tensor  and  consider  the  value  of  the  product  ak  mkl  m  when  ak 
lies  in  one  of  the  principal  directions  of  my.  Then  Ot*  m ki  =  a/,  where  Xm  is  the  eigenvalue  of 
mkl  corresponding  to  the  chosen  principal  direction.  The  product 

<*kmklZll  =  KalZll  (5-33> 

has  the  same  sign  as  the  eigenvalue  \m,  because  we  have  chosen  the  sign  of  Z]  in  such  a  way  that 
^  z  l  /  >  0.  (We  ignore  the  possibility  a;  z  j  /  =  0  because  it  arises  only  for  degenerate  elastic  behavior; 
note  that  ot/Zi/  =  + 1  for  isotropic  elasticity).  Because  mki  is  shear-like,  the  three  Xm  have  both 
positive  and  negative  eigenvalues  among  them,  so  that  the  product  a*  mk[  z j  /  (which  is  a  continuous 
function  of  ak)  must  pass  through  zero  as  a*  spans  all  possible  spatial  directions.  Therefore,  for  all 
practical  rate-independent  plasticity'  models,  there  exist  states  O  and  directions  ak  for  which  >  j 
(and/or  p\)  passes  through  zero,  (For  example,  in  the  isotropic  elastic  case  z1(  =  a,,  so  that  yt  and 

Y  p  y  p 

Pl  are  the  simple  quadratic  forms  -a^yay  and  -a,77i(yay,  respectively;  since  m,y  and  my  are 
shear-like,  these  forms  must  pass  through  zero  as  oc,  spans  the  directions  between  those  correspond- 

ing  to  the  maximum  and  minimum  principal  values  of  m,y  and  m,y).  It  is  clear  that  in  the 

non-associated  case,  m,y  *  m(y,  the  directions  corresponding  to  y\  =  0  will  be  different  from  those 
for  pi  =  0.  (This  can  be  rigorously  proven  by  means  of  a  straightforward,  but  lengthy,  argument). 
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As  a  consequence,  the  product p\y\ must  take  on  negative  values  (because  its  factors  pass 
through  zero  one  at  a  time).  Therefore  in  all  practical  non-associated  plasticity  models  there  will 
exist  states  O  and  directions  cr,  for  which  plastic  waves  can  propagate  with  a  sp'  *d  consistent 
with  M  >  A/f. 
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SECTION  6 

FAST  PLASTIC  WAVES  AND  NON-UNIQUENESS 


We  now  proceed  in  a  manner  similar  to  (Sandler  and  Rubin  1987)  in  order  to  show  that  the 

r 

existence  of  plastic  waves  with  M  >  M\  implies  non-uniqueness.  We  will  do  this  by  constructing 
multiple  (actually,  infinitely  many)  solutions  to  a  single  initial/boundary  value  problem.  Consider 
a  half  space  of  material  initially  at  rest  (at  yield  state  O)  but  subjected  to  a  prescribed  boundary- 
traction  which  varies  linearly  with  time  in  such  a  way  that  it  can  produce  a  response  involving  only 
an  elastic  (i.e.,  unloading)  primary  wave.  We  shall  show  that  this  problem  has  alternative  solutions 
(involving  plastic  response)  if  the  normal  a,  to  the  boundary  plane  lies  in  a  direction  for  which 
Pl  yi  <  0  in  state  O. 

The  general  nature  of  these  “spurious”  (but  nevertheless  mathematically  valid)  alternative 
solutions  is  shown  in  Figure  3.  In  this  figure  the  abscissa  is  a,  xt,  which  represents  the  position  of 
any  point  measured  from  the  boundary  plane  of  the  half-space,  a,  r,  =  0.  The  ordinate  is  time,  t, 

measured  from  some  arbitrarily  chosen  instant  t=0.  The  initial  conditions  at  r  =  0  are  o,;  =  cfy  and 

v,  =  0  with  the  boundary  traction  aj  o,y  at  the  equilibrium  value  a y  a°.  The  half-space  is  subjected 
to  the  traction  boundary  condition  (which  holds  along  the  ordinate  in  Figure  3), 
o.j  Gjj  =  tty  afj-Bpi  z\it  where  B  is  an  arbitrarily  prescribed  positive  constant.  As  shown  in  the 
figure,  the  alternative  solutions  that  we  are  about  to  construct  for  this  problem  consist  of  a  quiescent 
Zone  1 ,  followed  by  the  front  OF  (propagating  with  wave  speed  c)  of  a  unidirectional  plastic  loading 
wave  in  Zone  2,  followed  in  turn  by  a  transition  OV  (which  travels  with  speed  V)  to  elastic  behavior 
in  Zones  3,4,5  and  6.  The  lines  OP,  OS  and  OT  separating  the  elastic  zones  are  characteristics  whose 
slopes  are  the  inverses  of  the  three  elastic  wave  speeds.  For  simplicity,  we  will  construct  only 
continuous,  piecewise-linear  alternative  solutions,  i.e.,  those  for  which  the  velocity  vector  and  stress 
tensor  are  linear  functions  of  jc,  and  of  t  in  each  of  the  zones  shown  in  Figure  3. 

The  unidirectional  plastic  loading  wave  in  Zone  2  has  the  motion  direction  (5,  defined  by 
eqs.  (5.21,5.25).  For  a  piecewise  linear  solution  the  loading  rate  g"  in  Zone  2  is  a  positive  constant, 
so  that  g'  >  0  varies  linearly  with  distance  behind  the  wave  front.  From  eqs.  (5. 2,5. 5),  v,-,  and  otJ 
vary  linearly  with  distance  behind  the  wave  front, 

Vi  =  c$ig'=:czmiymg''s  (6.1) 
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Figure  6.1.  Characteristic  space-tune  construction  showing  spurious  waves  (Resulting  From 
Homogenous  Initial  /  Boundary  Conditions) : 

Zone  l  is  quiescent 
Zone  2  involves  plastic  loading 
Zones  3, 4, 5,  and  6  involve  elastic  unloading 
c  =  Plastic  Loading  Wave  Speed 

V  =  Speed  of  Elastic-Plastic  Interface 
c*  -  Primary  Elastic  Wave  Speed 
c\  =  Secondary  Elastic  Wave  Speed 
cj  =  Tertiary  Elastic  Wave  Speed 
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®ij  ~  ®ij  Ljjki  &k  Z-mt  y m  8  s 


(6.2) 


where  s  =  ct  -  ap  xp  is  the  distance  behind  the  wave  front.  Along  the  plastic -elastic  interface,  line 
OV  in  Figure  3,  the  velocity  vf  and  the  stress  a]j  may  be  obtained  by  setting  t  =  ap  xJV  to  get 

vf  =  c(~l)g"ymap  xp  zmi  (6'3) 

°ij  =  oy  ~  Ym  Lijki  ak  v-p  *P  zmt  (6,4) 


The  governing  equations  for  the  elastic  zones  consist  of  the  equation  of  motion, 

(6.5) 

dt  p  dxj 


and  the  constitutive  relation, 


^a‘j  _  r  1 


dvk  dvi 
dx[  dxf. 


dvi 

=  CW3Z 


(6.6) 


Let  Rm  and  Smbe  six  arbitrary  real  constants,  and  let  7^  =  Tp  be  another  six  real  constants  which 
satisfy  the  three  equations 


TijCLj^O 


It  is  easy  to  show  (by  direct  substitution)  that  the  functions 


vi  ~  (^m  f  Qp  *p  )Z-mi 


and 


®ij  =  +  ffij  "F  Rm  P Zmi  )  ®-pxp  +  Sm  C^i  CL^  Zm[  t 

identically  satisfy  the  governing  eqs.  (6.5, 6.6),  if  one  invokes  the  identity 

3  3(otp  xp )  3  _  d _ 

dxj  dxj  d(-xp  xp)  ~  i  d(ap  xp) 


(6.7) 


(6.8) 


(6.9) 


(6.10) 
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and  uses  eq.  (6.7).  Given  suitable  auxiliary  conditions  (such  as  initial  or  boundary  conditions  for 
v;  and  Oy),  the  relationships  (6.7,6.8,6.9)  form  a  system  of  twelve  equations  which  determine  the 
twelve  constants  Rm,  Smand7y.  We  can  eliminate  the  T y  by  multiplying  eq.  (6.9)  by  a;  and  using 
eqs.  (6.7)  and  (5.12,5.14)  to  get 

„  2,  p  (6.11) 

O,  (  Oy  -  Oy  )  =  }jRm  Pap  *p  +  Kti 

Now  eqs.  (6.8,6. 11)  form  a  system  of  six  equations  for  the  Rm  and  Sm,  with  the  7y  being  obtainable 
directly  from  eq.  (6.9)  after  Rm  and  Sm  are  determined.  The  new  system  can  be  further  simplified 
by  multiplying  each  equation  by  zRl-  (uncoupling  the  terms  corresponding  to  different  values  of  m ) 
to  get 

tRn  +  ot p  Xp  Sn  —  Vf  zni  (6. 12) 

ap  xp  Rn  +  (cf  )2  tSn  =  a.j  (Oy  -  Oy  )z„/p  (no  sum  on  n)  (6. 1 3) 

in  which  cE  s  ^/yV/f/p  is  the  n*  elastic  wave  speed. 

Along  line  OV  in  Figure  3,  where  CLpXp  =  Vt,  the  values  of  v,  and  Oy  in  eqs.  (6. 1 2,6. 13)  must 

equal  vt-  and  Oy  as  given  by  eqs.  (6.3, 6.4).  Therefore,  by  utilizing  eqs.  (5.1 1)  we  get  (after  some 
algebra). 


cv+(cE )2 

R»  =  -  — - 7r^c(c-V)g"yn  (no  sum  on  n) 


(6.14) 


"  ^-(4)2 


5„  = 


P  ~\p 


”  v’-fcj)2 


cg"y„  (no  sum  on  n) 


(6.15) 


in  Zone  3. 


In  the  remainder  of  the  elastic  region,  Rn  and  Sn  must  be  constant  except  along  the 
2  E  2 

characteristics  (ctp  xp  )  =  (c„  t)  ,  where  eqs.  (6.12,13)  become  singular.  Along  these  lines  the 
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theory  of  hyperbolic  partial  differential  equations  permits  the  “constants”  to  have  different  values 
in  the  two  regions  on  opposite  sides  of  a  characteristic  as  long  as  the  system  (6. 12,6. 13)  remains 

satisfied.  In  particular,  along  ap  xp  =  cn  t  the  system  becomes 

Rn  +  c%  Sn  =  v,  zn/t  (no  sum  on  n)  (6  1 6) 


Rn  +  cnSn~  0 ij  ( Cij  -  <5°  )zni/(pcEn  t)  (no  sum  on  n)  (6- 1 7) 

E 

Because  the  left  sides  of  eqs.  (6.16,6.17)  involve  only  the  combination  ( Rn  +  cn  S„),  the  individual 
values  of  the  constants  may  be  discontinuous  across  the  characteristic  line  if  the  combination  is 
preserved.  Let  Rn  and  Sn  denote  the  values  of  the  constants  behind  the  /1th  characteristic.  Then 


(6.18) 


V+Cn 


-c(c-V)g  y„  (no  sum  on  n) 


As  each  of  the  characteristics  OP,  OS  and  OT  in  Figure  3  is  crossed  in  passing  from  Zone  3  to  Zone 
6,  another  pair  Rn,  S„  jumps  to  Rn,  Sn.  For  Zone  6  eq.  (6.9)  implies 

aj  -G?j  )  =  (Rrn  Pap  Xp  +  Sm  )Zmi  for  0<apxp<  eft  (619) 

The  boundary  condition  for  the  problem  at  hand  is  a;  (ay  -  )  =  -Bpx  z\i  t  on  ap  xp  =  0 

for  all  t  >  0  (i.e.,  along  line  ON  in  Figure  3).  Because  the  zmi  are  linearly  independent  vectors,  eq. 
(6.19)  can  satisfy  this  boundary  condition  if,  and  only  if,  S2,  and  S3  are  zero  and 

SxM^-Bpx  (6.20) 


The  elastic  region  is  therefore  characterized  by  the  following  sets  of  constants: 


/?j,  SXt  /?2>  S2,  /?3  and  S3 

in  Zone  3 

R\ '  S\  /?2,  S2,  /?3  and  S3 

in  Zone  4 

/?l  Sj  1?2,  0,  ^?3  and  S3 

in  Zone  5 

/?l  S]  Rj,  0,  f?3  and  0 

in  Zone  6 
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At  this  point  we  have  completed  the  construction  of  a  set  of  alternative  solutions  to  the  simple 
initial/boundary  value  problem  posed  at  the  beginning  of  this  section.  However,  we  have  not  yet 
proven  that  the  stress  given  by  eq.  (6.9)  is  actually  consistent  with  the  assumed  elastic  behavior  in 
Zones  3, 4, 5  and  6.  In  order  to  do  this,  note  that  eq.  (6.9)  gives  the  rate  of  change  of  stress  in  elastic 
regions  as 


®ij  ~  Cijkl  ak  Z-ml 


(6.21) 


in  which 


Em  =  < 


cJ-  V2 


v2-<4)2‘*  Ym 


-Bp\/M\  for  m  =  1 


0 


for  m-  2  or  3 


} 


when 


when 


otp  Jtp  kXp  Xp 

,  *  s  t  s  _ 


(6.22) 


t> 


aP  xp 


-m 


(E, 


m 


a„Xp 

is  undefined  at  r  =  — KE~ ). 


-m 


The  loading  rate  relative  to  the  yield  surface,  /E,  is 

Y  .  y  Y 

nij  Gij  =  Em  fly  Cijtf  (X^  zmi  =  Em  rttij  tx^  zm[  —  ~Em  ym 

In  Zone  3  this  becomes. 


(6.23) 


4dij 


3 

Qacg "X 

m~\ 


c2-  V2 

v2-^2 


Pm  ym 
M-MEm 


(6.24) 


after  using  eqs.  (5.25,5.29).  The  sign  of  the  m  =  1  term  on  the  right  side  of  eq.  (6.24)  is  negative 
because  Q,  a ,  c  and  g"  are  all  positive,  c>V>cf>0  (so  that  the  expression  in  brackets  is  positive), 
and  pi  yj  is  negative.  Furthermore,  the  bracketed  expression  can  be  made  arbitrarily  large  in 

magnitude  for  m  =  1  by  choosing  the  value  of  V  to  be  only  slightly  larger  than  cf .  Because  cf  and 
C3  are  less  than  c\ ,  the  m  =  1  term  will  then  dominate  eq.  (6.24)  and  determine  its  sign.  Therefore 
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values  of  V can  be  chosen  for  which  the  loading  rate  is  negative  in  Zone  3  (and  therefore  consistent 
with  the  presumed  elastic  behavior  of  Figure  3). 


In  Zones  4,  5  and  6  the  loading  rate  is,  according  to  eqs.  (6.23,6.22), 

nJ<jy  =  5pi>i/A/f  +  iVg" 


(6.25) 


in  which 


n  V  Pm  y'm 

,  £\2  M  mE 
m=2  *  ~  (cm) 


Qac 


c2  -  V2  Pi  >3 
V2  -  (cf )2  M  — A/f 


in  Zone  4 


in  Zone  5 


0 


in  Zone  6 


(6.26) 


Now,  the  value  of  g"  can  be  chosen  small  enough  so  that  the  first  term  on  the  right  side  of  eq.  (6.25) 
dominates  the  second.  Then,  because  B  and  A/f  are  positive  and  p\  y\  is  negative,  the  loading  rates 
in  Zones  4,  5  and  6  are  negative  (and  the  response  is  elastic  in  those  zones  as  required  in  Figure  3). 

Summarizing,  we  have  constructed  (and  confirmed  the  validity  of)  multiple  solutions  to  a 
simple  initial/boundary  value  problem.  These  solutions  represent  alternatives  to  the  simple  solution 
(involving  only  primary  elastic  unloading  waves  and  corresponding  to  g”  =  0).  The  existence  of 
such  multiple  solutions  (involving  appropriate  values  of  V  and  g"  >  0)  proves  that  non-associated 
flow  leads  to  non-uniqueness.  Because  the  plasticity  formulation  analyzed  here  is  completely 
general,  it  extends  over  the  entire  class  of  rate-independent  plasticity  models  used  for  practical 
ground  shock  calculations. 
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SECTION  7 
CONCLUSION 


Many  ground  shock  calculations  performed  for  the  Defense  Nuclear  Agency  utilize  rate-in¬ 
dependent  plasticity  models  with  non-associated  (i.e.,  non-normal)  flow  rules  to  represent  the 
behavior  of  geological  materials.  In  such  models  the  yield  surface  is  used  to  limit  the  shear  and/or 
tensile  stresses,  while  the  flow  rule  is  used  to  independently  define  the  direction  of  plastic  straining 
(thereby  limiting  the  dilatancy).  The  behavior  of  constitutive  formulations  of  this  kind  in  dynamic 
applications  has  been  analyzed  in  this  report.  By  means  of  a  “constructive”  proof,  it  has  been 
demonstrated  mathematically  that  such  models  lead  to  multiple,  and  therefore  spurious,  solutions 
in  wave  propagation  situations;  in  other  words,  uniqueness  of  solution  breaks  down  for  dynamic 
initial/boundary  value  problems. 

Although  a  simple  class  of  multiple  solutions  was  presented  to  complete  the  proof  for  a  very 
simple  type  of  problem,  the  conclusion  is  quite  general;  the  proof  makes  no  specific  modeling 
assumptions,  except  rate-independence,  shear-like  plasticity  and  non-singular  elastic  behavior  (and 
is  therefore  relevant  to  all  of  the  ground  shock,  and  some  of  the  structural,  plasticity  models  used 
in  practice).  The  loss  of  uniqueness  is  found  to  be  closely  related  to  the  energy  generating  properties 
of  the  models  according  to  D’iushin’s  postulate  (even  though  the  second  law  of  thermodynamics 
may  be  satisfied).  The  result  implies  that  normality  of  flow  is  a  necessary  (as  well  as  a  sufficient) 
condition  for  uniqueness  in  dynamic  applications  of  any  rate-independent  plasticity  model. 

Obviously,  these  findings  imply  an  inherent  lack  of  robustness  in  numerical  analyses  based 
on  non-associated  plasticity.  Because  the  use  of  these  models  is  widespread  in  ground  shock  (and 
structural)  calculations,  the  fact  that  they  can  produce  multiple  (and  therefore  spurious)  solutions  is 
significant;  it  casts  serious  doubt  on  the  validity  of  the  specific  results  of  any  (and  every)  calculation 
based  on  this  kind  of  constitutive  representation.  Therefore  it  is  strongly  recommended  that  DNA 
promptly  abandon  its  reliance  on  such  models  in  all  work  performed  for  the  Agency.  To  describe 
the  situation  most  bluntly,  these  models  do  not  provide  a  rational  basis  for  the  computations  needed 
for  prediction,  design  and/or  analysis  applications.  In  all  cases  for  which  an  associated  flow  rule  is 
deemed  inadequate  to  fit  observed  material  behavior,  an  approach  other  than  rate-independent 
plasticity  must  be  sought  in  order  to  represent  such  behavior  in  a  reliable,  self-consistent  and  rational 
manner. 
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In  view  of  the  fact  that  non-associated  plasticity  models  have  already  been  extensively  used 
in  analyses  for  DNA,  it  is  recommended  that  some  effort  be  made  to  quantitatively  assess  the 
consequences  of  their  use.  In  particular,  the  extent  to  which  the  model  assumptions  and  parameters 
can  affect  the  numerical  solutions  to  various  types  of  problems  should  be  investigated. 
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APPENDIX 

DERIVATION  OF  THE  STRAIN  INCREMENTS 
SATISFYING  INEQUALITY  (4.7) 


In  this  Appendix  we  present  the  entire  six-dimensional  set  of  strain  increments  which 
satisfy  eq.  (4.7)  whenever  the  flow  rule  is  non-associated.  Consider  the  ratio  of  each  of  the 
components  ny  to  the  corresponding  component  ny.  Either  these  ratios  are  all  the  same  (i.e. 
independent  of  i  and  J)  or  two  (or  more)  of  them  differ.  If  the  former  possibility  holds,  we  can  let  r 
be  the  common  ratio,  so  that  ny  =  r  rty  for  all  ij.  But  in  that  case  ny  riy  =  r  ny  ny,  so  that  r  =  ±1 . 

When  r  =  +1  the  flow  rule  is  associated,  so  this  situation  is  not  relevant  here.  For  r  =  -I,  ny  =  -ny, 
and  eq.  (4.7)  holds  for  any  A £*/  satisfying  eq.  (3.6).  The  only  remaining  possibility  in  the  case  of 

non-associated  flow  is  that  two  (or  more)  of  the  six  ratios  ny  /ny  are  different.  Then  we  can  select 
at  least  two  pairs  of  subscripts  m,n  and  p,q  in  such  a  way  that 

nmn  (A- 1 ) 

— r~  *  (no  “sum”  on  m,n  or  on  p,q  ) 

nmn  npq 

Let  h\,  h2,  h2  and  h4  be  four  arbitrary  real  numbers.  Form  a  symmetric  tensor  hy  which  has 
zero  components  in  the  m,n  and p,q  positions  and  has  h\,  h2,  h2  and  h4  as  the  other  four  independent 

Y  Y  P  P 

components  (in  any  order).  Also  define  h  &  nyhy,  h  =  nyhy  and  D  as 

D  =  nlmKpq  -  nYpqnPmn*  0  (A.2) 

Then  two  additional  arbitrary  real  numbers,  p  and  v,  together  with  h\,  h2,  h2  and  h4  define  a 
six-parameter  family  of  “elastic  stress  increments” 

Cykl^-kl s  M-[(  l-v)«y~(l+v)ny]+/iy  +(h  n^-h  nffm)Uypq—(h  npg—h  npq)Uymn  (A.3) 


in  which  Uijmn  *  ~  (6 >im6jn  +  8^5^).  By  multiplying  eq.  (A.3)  by  ny  and  ny,  respectively,  it  is 


straightforward  to  verify  that  eq.  (4.7)  is  satisfied  for  all  p  >  0  and  Ivl  < 


y  p 

1  -  nq  nq 

Y  P‘ 
1  +  ny  ny 
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